Quantum state diffusion, localization and computation 



Riidiger Schack, Todd A Brim and Ian C Percival 
Department of Physics 
Queen Mary and Westfield College, University of London 
Mile End Road, London El 4NS, England. 

February 1, 2008 



Abstract 

Numerical simulation of individual open quantum systems has proven advantages 
over density operator computations. Quantum state diffusion with a moving basis 
(MQSD) provides a practical numerical simulation method which takes full advantage 
of the localization of quantum states into wave packets occupying small regions of 
classical phase space. Following and extending the original proposal of Percival, Alber 
and Steimle, we show that MQSD can provide a further gain over ordinary QSD and 
other quantum trajectory methods of many orders of magnitude in computational 
space and time. Because of these gains, it is even possible to calculate an open quan- 
tum system trajectory when the corresponding isolated system is intractable. MQSD 
is particularly advantageous where classical or semiclassical dynamics provides an ade- 
quate qualitative picture but is numerically inaccurate because of significant quantum 
effects. The principles are illustrated by computations for the quantum Duffing oscil- 
lator and for second harmonic generation in quantum optics. Potential applications 
in atomic and molecular dynamics, quantum circuits and quantum computation are 
suggested. 
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1 Introduction 



Most quantum systems are not even approximatedly isolated, but open, so that they are 
signficantly affected by the environment. This interaction is important for atoms and 
molecules in gaseous or condensed matter environments, which broaden spectral lines. It 
affects the motion of molecules and the rates of chemical reactions. It is important for 
signalling near the quantum limit, where the environment produces the noise through which 
the signal must be detected. And it is important in quantum optics, where it produces the 
dissipation that destroys coherence. 

Quantum state diffusion with a moving basis (MQSD) is a method of representing and 
computing the evolution of individual open quantum systems. It has already been used 
by Percival et al. |39| to analyse the motion of a particle in a Penning trap. Here we 
provide a general theory of the method, and provide a guide as to when it should be used 
in preference to other methods. Second harmonic generation in optics and the quantum 
Duffing oscillator are used as illustrations. 

Because stochastic environmental fluctuations affect the evolution of an individual open 
quantum system, it is represented traditionally by a density operator p, which satisfies 
a linear master equation. No attempt is made to represent the evolution of individual 
pure states explicitly. This approach is adequate when the master equation has analytic 
solutions, or when the number of basis states N of Hilbert space required for numerical 
solution is not too large. But for large iV it often breaks down in practice long before 
the corresponding numerical solution of the Schrodinger equation, because the number of 
elements of a density matrix increases as N 2 . 

Numerical simulation of individual open quantum systems, represented by pure states 
which move along quantum trajectories, has proven advantages over density operator com- 
putations. Quantum state diffusion (QSD) provides such a numerical simulation, in which 
each state diffuses continuously in the state space and satisfies a nonlinear Langevin-Ito 



diffusion equation, determined uniquely by the master equation as shown in [25] and de- 
scribed in section |2|. This diffusion often produces a localization of quantum states into 
wave packets that occupy small moving regions R of classical phase space. Such localization 
is a special characteristic of QSD that is usually absent in other quantum state simulation 



methods [16]. 

For a system of / freedoms, a Planck cell of volume {2^h)^ corresponds to one quantum 
state. The number of basis states used to represent a system need not be much greater than 
the number of Planck cells in the region R, provided that the basis follows the motion of the 
wave packet in phase space. Quantum state diffusion with a moving basis (MQSD) provides 
a practical numerical simulation method which takes full advantage of the localization, by 
referring the quantum state to a moving origin (q,p) in phase space. This origin lies at the 
phase space centroid of the quantum state, determined by the current quantum expectations 
(Q) and (P). 

Numerical methods for solving the time-dependent Schrodinger equation for an isolated 
system using moving wave packets have long been used in chemical physics [[31], |(| |32[ . 
But as the ordinary Schrodinger equation disperses wave packets instead of localizing them 
in phase space, the applicability of these wave packet approaches is very restricted. 
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For some systems, computation using MQSD is orders of magnitude more economical 
in computer storage space and computation time than other quantum state simulation 
methods, or the solution of master equations for the density operator. We give examples in 
which it is very difficult to see how any other current method of numerical solution could 
be used. 

Section ^| presents the basic QSD equations and their derivation. The problems of 
the choice of boundary between system and environment are described. There is a brief 
comparison with other quantum trajectory methods, sometimes called quantum jump or 
relative state methods || |7], [14], T5|. Localization is defined and discussed in section |3|, 



with reference to localization theorems and numerical examples. This is followed by the 
definition of the moving basis, using excited coherent states, and the derivation of the 
MQSD equations. 

Section [| applies MQSD to two challenging examples, the Duffing oscillator and second 
harmonic generation. These two examples show how effective MQSD can be in bridging 
the gap between those quantum problems where the number of basis states required in a 
fixed basis is relatively small, and the quasiclassical limit where the number of fixed basis 
states would be so large as to rule out any practical use. This section is completed by a 
crude analysis of comparative computing times for the solution of an open system problem 
using MQSD and the numerical solution of the Schrodinger equation for a similar isolated 
system. 

Section [5] concludes with a comparison of methods for open systems and some recom- 
mendations for their use, followed by the prospects for using MQSD in various applications. 



2 Quantum state diffusion (QSD) 

Quantum state diffusion represents the evolution of a quantum system through a corre- 
spondence between the solutions of the master equation for the ensemble density operator 
p and the solutions of a Langevin-Ito diffusion equation for the normalized pure state vector 
of an individual system of the ensemble [f25|j . 

An analogy is helpful. A solution of the Langevin-Ito equation for the motion of an 
individual Brownian particle in position space represents a single member of an ensemble 
whose distribution function satisfies the corresponding Fokker-Planck equation. Similarly, 
a solution of the QSD equation for the diffusion of a pure quantum state in state space rep- 
resents a single member of an ensemble whose density operator satisfies the corresponding 
master equation. In each case the effect of the environment can be represented either by the 
stochastic evolution of an individual system or by the deterministic evolution of the distri- 
bution. For Brownian motion the evolution of the individual system gives a more detailed 
picture of what happens to an individual particle than the evolution of the distribution 
function. Similarly, for open quantum systems, the evolution of the individual pure states 
of QSD gives a more detailed picture of what happens to an individual quantum system 
than the evolution of the density operator. This is particularly important for applications 
like single particle traps, quantum noise in gravitational wave detection, quantum circuits 
and quantum computers. 
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In QSD, quantum expectations (. . .) for individual systems, and ensemble means M are 
distinct. The traditional quantum expectation TrpG of an operator G for a mixed state of 
an open system is equivalent in QSD to an ensemble mean over the quantum expectations 
of the pure states \ip). That is 

TrpG = M(G) = M(i/)\G\ip). (1) 
If the master equation has the standard Lindblad pi form 



9 = ~^[ H ' 9] + ( L mP L m - |L^L m p - ipL|„L m ) , (2) 



Til 



then the corresponding QSD equation is a nonlinear stochastic differential equation for the 
normalized state vector of the ensemble, whose general differerential form is 

\dll>) = ~im)dt + Em «l4>Lm " kMn^m ~ § (LjJ <L m » \^)dt 

+ E m (L m - (L m ))\if>)d£ m , (3) 

where H is a Hamiltonian and L m are Lindblad operators which represent the effect of the 
environment on the system in a Markov approximation. 

The first sum in (||) represents the 'drift' of the state vector in the state space and 
the second sum the random fluctuations. The d^ m are independent complex differential 
random variables, with normalized independent white noise in their real and imaginary 
parts, leading to an isotropic Brownian motion or Wiener process in the complex £ m -plane. 
These satisfy the conditions 

Md£ m = 

Md£ n d£ m = 0, MdCdU = S nm dt, (4) 

where M represents a mean over the ensemble. The complex Wiener process is normalized 
to dt, so the independent real and imaginary Wiener processes are each normalized to di/2. 



This is the same normalization as in 26, 27, 381, but is different from that of [Eql. The 



distribution (4) is invariant under unitary transformations in the linear space of the d£ m . 

(L m ) = (V'|L m |'0) is the quantum expectation of L m for state \ip). The density operator 
is given by the mean over the projectors onto the quantum states of the ensemble: 

P = M|V><# (5) 

It can be verified that if the pure states of the ensemble satisfy the QSD equation (^|), 
then the density operator @ satisfies the master equation (D). There are many diffusion 
equations for pure states which give the same master equation for the density operator. 
The uniqueness of the QSD equations follows from a principle of unitary invariance in 
operator space. 

This is most clearly illustrated by a QSD equation with a single Lindblad operator L. 
In that case unitary transformation in operator space is multiplication by a scalar phase 
factor u of modulus unity. Then it is obvious by inspection that if L is replaced by uL, 
the master equation, and hence the solution of the master equation, is unchanged. In the 
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corresponding QSD equation, the replacement of L by uL is the same as the replacement 
of d£ by wd£, but since the distribution of the elementary differential fluctuations d£ is 
invariant for multiplication by a phase factor, the QSD equations are also unchanged. This 
is not true if the fluctuations are real or otherwise do not satisfy unitary invariance in 
fluctuation space, as in [T9], 122]. The QSD equations are the only diffusion equations that 



respect this unitary invariance property in the one-dimensional operator space, and they 
are in this respect unique up to a physically irrelevant external time-dependent phase factor 
for the state vector. 

This result generalizes to an arbitrary number of Lindblad operators if u is taken to be 
a general unitary transformation in the linear space of the Lindblad operators. The general 



result was first given for QSD equations in p5 |, following the statement of an invariance 



principle by Diosi B, and the derivation for a Fokker-Planck equation in state space in 



In either the traditional or QSD formulation for open systems, there is always a problem 
as to where to put the boundary between system and environment. If the system is made 
too small, then important effects are neglected, and errors are made. If it is made too 
big by including too much of the environment, then the dimension N of the basis state 
space becomes too large for the equations to be solved. In practice there is a compromise. 
Often, in addition, the problem is simplified by approximating the effects of complicated 
environments by simple operators. 

Because of the localization property discussed in the next section, QSD has the property 
that it has no need of a separate measurement hypothesis, as shown in detail with examples 
in [2J|, |26], p7| . Measuring apparatus is just one type of environment, and its effects can be 
represented by simple operators. Alternatively, measurement can be represented in detail by 
pushing the boundary between system and environment out so far that the system includes 
the measuring apparatus. This is useful in establishing the representation of measurement 
by operators, but results in a system too complicated to solve directly. Localization and the 
representation of measuring apparatus are needed for comparison with the relative state or 
quantum jump methods. 

QSD theory followed from research that was motivated by the desire to find an explicit 
physical representation of the measurement process. Following pioneering work of Bohm 
and Bub |3J and Pearle [^, |3E], Gisin |2D| introduced a simple example of quantum state 



diffusion with real fluctuations that was generalized by Diosi and Gisin [pi]] . The complex 



Ito form of QSD was introduced in [p5| . The detailed QSD theory and its applications are 
described in |25], ^7], |38fl . Diffusion in the space of quantum states also appears in 
connection with the theory of continuous quantum measurements as shown, for example, 
in the many references in ]!], Q . 

Gisin and Percival pl| used QSD to describe a quantum jump experiment. Goetsch and 
Graham [^3 used it to describe some nonlinear optical processes. Garraway and Knight O] 



compared QSD and quantum jump simuations for two-photon processes, and in |18j they 
compared the phase space picture of QSD and quantum jumps, showing that the former 
gives localization and the latter does not. Spiller and his coworkers applied QSD to thermal 



equilibrium |44j], studied chaos in a simple open quantum system E| and investigated open 



angular systems, such as quantum capacitors and rotors ffEfl. Gisin |2j| investigated the 
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Heisenberg picture for QSD. Further examples are given in |25|, |26|, 



G71 



EDI 



3 Localization and the moving basis 

For a quantum system in a pure state, localization refers to dynamical variables and oper- 
ators like position, momentum, energy and angular momentum which have classical equiv- 
alents. It is helpful to picture the localization in classical phase space. There are many 
phase space representations of quantum systems, such as the Wigner distribution, but these 
are not necessary to obtain a classical picture of the phase space localization of a quantum 
system. For that, it is only necessary to consider the expectations and variances of quan- 
tum dynamical variables, and then to picture these expectations and variances as if 'they 
represented the expectations and variances of classical quantities. 



The general theory of localization is treated in |25] and in |38| . The Schrodinger evo- 
lution of a system usually produces (^localization or dispersion. The interaction of the 
system with its environment, by contrast, produces localization. 

It is convenient to consider these competing processes separately before considering 
them together. The dispersion is well-known, since it occurs in isolated systems. The 
opposite extreme is an open system which interacts with environment so strongly that the 
Schrodinger evolution can be neglected. This is a wide open system, and the theory of 



localization for these systems has been treated in detail in J38 



For simple systems, earlier papers E3, E7fl had built up a picture in which interaction 



with the environment produced localization, sometimes to an eigenstate of an operator 
corresponding to a surface in phase space, but more commonly to a state which is localized 
to a wave packet in phase space whose Heisenberg indeterminacy products are of the order 
of h. In [B3| that picture was confirmed and extended. A general theory was presented, 
lower bounds were put on rates of self-localization, and bounds were put on asymptotic 
states. 

For a pure state \if)), the expectation of a selfadjoint operator G is denoted 

(G) = MG|V>> (6) 

and the variance of the operator G is 

a 2 (G) = (G 2 ) - (G) 2 . (7) 

The localization A is defined as the inverse of the mean of the variance for the ensemble: 

A = (Mcr 2 (G)) -1 . (8) 

The simplest theorem is for a wide open system in which there is only one self-adjoint 
Lindblad operator L in the state diffusion equation. It is shown that the rate of localization 
of L is at least as fast as 2 (1331, section Ef): 



dA/dt > 2. (9) 
Because Lindblad operators have dimension time -5 , the rate is dimensionless. 
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This localization is towards a surface in phase space defined by the dynamical variable 
L, and is characteristic of interaction with apparatus that measures L. More general and 
more interesting is the case where there are at least 2f operators, where / is the number of 
freedoms, and the surfaces for the corresponding dynamical variables define points in phase 
space. In that case at least / pairs of operators will not commute. The simplest example 
is a one freedom system with two conjugate variables X, Y whose operators satisfy 

pC,Y] = ift. (10) 

In this case it is shown that the state localizes and asymptotically approaches a wave 
packet with minimum Heisenberg indeterminacy product. It is also shown that operators 
that are not selfadjoint, such as annihilation and creation operators, can localize to wave 
packets. In the more realistic examples that we study numerically in the following section, 
the localization has to compete with the dispersion due to the Schrodinger evolution, the 
wave packet is more dispersed, and its dispersion often varies considerably as a function of 
time. 

For general open systems there is a physical competition between the dispersion due to 
the Hamiltonian and the localization due to the state diffusion. There is no general theory 



for this, but there are many numerical examples of localization, such as [25| for the forced 



damped oscillator and |27| for localization in one well of a double well, and phase space 



localization by position localization and Hamiltonian coupling of position and momentum. 
Halliwell and Zoupas have provided a general theory for the evolution of Gaussian wave 
packets for an important model, generalizing a result of Diosi ||. 

From the theorems and the numerical examples it would appear that the localization 
of wave packets to relatively small regions of phase space is the norm, so that the limiting 
behaviour on a classical scale is the direct representation of classical states as points in 
phase space, rather than the more abstract surfaces defined by action functions that satisfy 
the Hamilton- Jacobi equation. 

For our purposes the most important consequence of the localization of quantum trajec- 
tories around phase space trajectories is, that by continually changing the basis, it is often 
possible to reduce the number of basis states needed to represent the wave packet by many 
orders of magnitude. If a wave packet is localized about a point (q,p) in phase space far 
from the origin, it requires a great many of the usual number states \n) to represent it. But 
fewer excited coherent basis states \q,p, n) = ~D(q,p)\n), are needed, with corresponding 
savings in computer storage space and computation time. These states are defined using 
the coherent state displacement operator, 

V(q,p) = ex V UpQ-qP). (11) 



h 

The separation of the representation into a classical part (q,p) and a quantum part \q,p, n) 



is called the moving basis, or, as in [39], the mixed representation 



In this basis, the usual creation and annihilation operators are modified: 

a = a.(q,p) + (q + zp)I/V2, a f = eJ(q,p) + (q - ip)l/y/2, (12) 
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where a(q, p) is the local annihilation operator (with (q, p) as the origin) and I is the identity 
operator. The effect of the local operators on the excited coherent states is given by: 

a(g,p)|g,P,n) = Vn\q,p,n- 1), a f (q,p)\q,p,n) = Vn + 1 \q,p,n + l). (13) 
Similarly, 

Q = Q(q,p) + qI, P = P(q,p)+pI, (14) 

where 

Q(q,p) = (a(g,p) + a t (g,p))/v / 2, P(g,p) = -i(a(g,p) - a!(q,p))/V2. (15) 

Care must be taken to avoid ambiguity in the external phase factors. Normally, different 
displacement operators do not commute; nor is the product of two displacement operators 
a standard displacement operator. In both cases there is an additional phase factor: 

B(q',p')\q,p,n) = \q + q', p + p f , n)e^'~^ 2 . (16) 

In order to retain an unambiguous relation to the standard fixed basis states, this additional 
phase factor must be removed. Fortunately this is simple, as it the same for all the moving 
basis states. 

The numerical algorithm follows directly. As the integration proceeds, the phase point 
expectation 

(6q,5p) = ({Q(q,p)),(P(q,p))) (17) 
drifts away from zero. The basis is then shifted, 

\^)^-D(-Sq,-SpM, (18) 

and the classical phase point (q,p) adjusted: 

(q,p) -> (q + Sq,p + 5p). (19) 

The time between basis adjustments can be optimized in various ways, depending on the 
problem, the degree of localization, and the difficulty of carrying out the basis change. 
In principle there is another way to make the basis change. If the phase space point 

(g,p) = «Q>,(P» (20) 

and state \ip) relative to (q,p) as origin are used to define the state, then we can include 
the change of basis in the evolution equation for This leads to a set of simultaneous 
differential equations in q, p, and in integrating these equations, both the QSD evo- 
lution and the basis shift are done automatically. This simultaneous moving basis method 
is attractive for a number of reasons; the states always remain in an optimally 'localized' 
basis, and therefore can take the best advantage of the QSD localization effects to mini- 
mize the number of necessary basis states. There is no need for an additional step for the 
basis change. In practice, the simultaneous moving basis method results in a considerable 
increase in programming complexity, and may have numerical stability problems. As yet 
the method has not proved itself sufficiently, so we have not used it. 
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Implementing the moving basis algorithm on a computer is straightforward. Each set of 
normalized ffucutations d£ m determines a quantum trajectory through Eq. (]3|), which can be 
simulated using discrete time steps St, using a Runge-Kutta algorithm for the deterministic 
part, and an Euler method for the stochastic part (but see also [47]). Suppose that at time 
t — to the state \ip(t )) is represented in the basis \qo,Po,n), centered at 

(?o,po) = (mo)\Qmt )), (v>(to)iPiv(to)». (21) 

Then after one discrete time step, the expectations in this basis shift to 

(q' ,p' ) = «^(to + 6t)\Qm + St)), (iP(t + St)\P\i>(t + St))) + (q ,p ). (22) 

The computational advantage of a small number of basis states is then retained by 
changing the representation to the shifted basis \qi,Pi,n) centered at q x and p x . This shift 
in the origin of the basis represents the elementary single step of the moving basis of MQSD. 

The components of \ip(t + St)) can be computed using the expressions given above. The 
computing time needed for the basis shift is of the same order of magnitude as for computing 
a single discrete timestep of Eq. (Q). Shifting the basis once every discrete timestep could 
therefore double the computing time, depending on the complexity of the Hamiltonian and 
the number of degrees of freedom. On the other hand, the reduced number of basis vectors 
needed to represent states in the moving basis can lead to savings far bigger than a factor 
of 2. 

In one example of second harmonic generation described in the following section, two 
modes of the electromagnetic field interact. Using the moving basis reduces the number of 
basis vectors needed by a factor of 100 in each mode. The total number of basis vectors 
needed is thus reduced by a factor of 10000, leading to reduction in computing time by a 
factor of 10000/2 = 5000. Futhermore, the fixed basis would exceed the memory capacity 
of most existing computers. 

The QSD equation (|3]) can contain both localizing and delocalizing terms. Nonlinear 
terms in the Hamiltonian tend to spread the wave function in phase space, whereas the 
Lindblad terms localize. Accordingly, the width of the wave packets varies along a typical 
trajectory. We use this to reduce the computing time even further by dynamically adjusting 
the number of basis vectors. Our criterion for this adjustment depends on parameters e < 1, 
the cutoff probability, and -/V pa d, the pad size, which represents the number of boundary basis 
states that are checked for significant probability. We require the total probability of the 
top A/p^ states to be no greater than e, increasing and decreasing the number of states 
actually used accordingly, as the integration procedes along the quantum trajectory. 



4 Examples: The Duffing oscillator and second har- 
monic generation 

The quantum mechanics of systems whose classical limit exhibits dissipative chaos is an 
interesting problem. Dissipation is relatively difficult to treat in quantum mechanics. The 
best of the commonly used techniques is the solution of the master equation, but as we have 
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pointed out, solving the master equation numerically can be an extremely difficult problem. 
Perhaps because of this, quantum dissipative chaos has been neglected by comparison 
with quantum hamiltonian chaos. Clearly the QSD method for open systems is eminently 
suitable for such dissipative systems. Recently Spiller and Ralph |45| applied QSD to a 
dissipative chaotic system. 

The Duffing oscillator is particularly appropriate for applying MQSD, as the classical 
system has been widely studied |29| . This system has also been treated quantum mechani- 



cally in the decoherent histories formalism B 0j, which has been shown to have connections 



to QSD ||10|| . This oscillator consists of a particle moving in one dimension in the two- welled 
potential 



4 2 
X X 



V(x) = - - -. (23) 

Dissipative chaos can be produced by adding both dissipation (of the form — 2Tx) and a 
periodic driving force (of the form gcos(t)). 

By using either the master equation or QSD, it is relatively simple to calculate the 
evolution of this system far from the classical limit, for example when tl — 1. But it is 
more important to study the classical limit of quantum chaos, to see how the relatively 
well understood properties of classical chaos appear in quantum systems. This limit is 
conveniently represented by decreasing h, which increases the number of fixed basis states 
required. With h — 1, about 10 states are needed to simulate the system with accuracy. 
On approaching the classical limit, the number of states quickly becomes impractical. In a 
semiclassical regime with h « 10 -4 , more than 10,000 states are needed; the density oper- 
ator method would require storage and computation with a prohibitive 10 8 real numbers. 
This, like most nonlinear problems, is exacerbated by the fact that the potential needs 
much more computation than, for example, the simple harmonic oscillator. 

Using MQSD, however, this problem becomes tractable. In a chaotic system, the de- 
localizing forces are particularly strong in the neighbourhood of hyperbolic fixed points, 
where the dynamics spreads the wave packet in phase space. However, the localizing effects 
of the environment always predominate. MQSD picks out a good local time-dependent ba- 
sis, and the method rarely requires more than twenty states, usually about ten, just as in 
the small scale quantum limit. This is easy to see by examining figure |I|. 

Our second example is frequency doubling or second harmonic generation, which is a 
standard process in quantum optics. The system consists of two optical modes of frequency 
uji and u 2 — 2uj\ which interact in a cavity driven by a coherent external field with frequency 
ujf ~ l)\ and amplitude /. The cavity modes are slightly damped and detuned, with 
detuning parameters 8\ = u>i — ujf and <5 2 = 0J2 — 2a;/. 

The Hamiltonian in the interaction picture is 

H = fttf^ai + M 2 a+a 2 + ihf{a{ - ai) + ^(a^a, - a?aj) , (24) 

where &i and a 2 are the annihilation operators of the two cavity modes, and x describes 
the strength of the nonlinear interaction between them. Damping of the two cavity modes 
is described by the Lindblad operators L x = y / 2l^"a 1 and L 2 = y / 2K 2 "a 2 . The factors of \[2 
are a consequence of normalization conventions used in the master equation Eq. (0) which 
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differ from those commonly used in quantum optics. The master equation for this problem 
first appeared in |12], |j~3f . 

A direct numerical solution of the master equation for this problem is difficult because, 
in the Fock basis, the dimension of the effective Hilbert space is equal to the product of the 
photon number cutoffs in both modes and easily becomes very large. The problem becomes 
intractable even for moderate numbers of photons in each mode. Earlier treatments of this 



problem include |fLl|, fLy, [43|, [2£ 

The QSD method was employed for studying second harmonic generation by Gisin and 
Percival ]25| and Goetsch and Graham [ 2q| . Both used a fixed Fock-state basis, and rather 



limited photon numbers. Gisin and Percival improved the method slightly by using 
a lower as well as an upper cutoff. More recently, Zheng and Savage |^8| applied the 
quantum jump method to the chaotic regime of second-harmonic generation for large 
photon numbers. Using hundreds of hours on a 32-processor supercomputer, they succeed 
in simulating cases where the expectation of the photon number is of the order of 200 in 
each mode. Their Fock basis has 512 basis states in each mode, i.e., a total of 250000 basis 
states, a tour-de-force near the limits of their method. 

To illustrate the power of MQSD, we have computed trajectories for second harmonic 
generation with parameters similar to the ones used by Zheng and Savage ^8|, but leading 



to photon numbers on average six times as large. For this problem, Zheng and Savage would 
have needed approximately 36 x 250000 ~ 10 7 basis states. Using the moving basis with 
an adaptive basis size, we computed a single trajectory in a few hours on a PC. The results 
shown in Fig. || were obtained using a cutoff probability e = 10~ 3 . The total number of 
basis states needed in the simulation varied between approximately 50 and 1000. Changing 
the cutoff probability to e = 10~ 2 lowered the number of basis states needed to between 
approximately 40 and 200, clearly a strong dependence on e. Surprisingly, the choice of 
cutoff probability affected the simulated expectation values only very slightly. 

It is interesting to compare estimates of computation space and time of a single MQSD 
run for an open sytems and the numerical solution of the Schrodinger equation for the 
same system isolated from its environment. The former has diffusion and drift terms that 
increase the computation time for a single step of integration, but the localization reduces 
significantly the size of the basis needed for computation. The latter has fewer terms in 
the equation, but the Schrodinger dispersion increases the number of states needed for the 
representation. Clearly, for systems with more than one degree of freedom the advantage of 
MQSD is even greater, the savings in number of states going like the power of the number 
of freedoms. We can make this semi- quantitative by the following analysis: 

Consider computing times. For each computation, let Nj be the arithmetic mean of 
the number of basis states required for the jth freedom. Let N be the geometric mean 
of all the Nj, which is is an appropriate measure for the number of basis states for one 
freedom. Let / be the number of freedoms. Then an estimate of the computing time T comp 
is given by KN*, where the constant K is much larger for MQSD but N is much larger for 
Schrodinger. The ratio of computing times is then 

T comp (Sch) _ K(Sch) r iV(Sch) ]f 



T comp (MQSD) if (MQSD) LjV(MQSD) 
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The values of these numbers are highly problem- dependent. Clearly, for a problem with 
a complicated Hamiltonian and one simple Lindblad operator (as in the Duffing oscillator), 
the difference in time between the QSD and Schrodinger calculations will be much less 
than it would be for a simple Hamiltonian with several Lindblad operators (as in second 
harmonic generation). 

One can estimate the values of these numbers for the problems described. For the Duff- 
ing Oscillator, K (Sch) / K (MQSD) « 2/3, iV(Sch)/iV(MQSD) « 1500, and / = 1. From 
(|25| ) we expect that the ratio of computing time between MQSD and the Schrodinger equa- 
tion is roughly 10~ 3 . In the case of second harmonic generation, we have K (Sch) / i^(MQSD) m 
1/2, iy(Sch)/iV(MQSD) ps 100 and / = 2, yielding a ratio of computing times of approxi- 
mately 2 x 10~ 4 . Clearly, the advantage for multi-freedom problems is very great. 



5 Conclusions and prospects 



For the simplest open system problems, where analytic solutions are available, or the num- 
ber N of basis states that are needed for the solution of the master equation is small, a 
density operator method is to be preferred. This method may be the only one avaibable 
when very high accuracy is needed, because of the inevitable statistical errors in Monte 
Carlo simulations. 

Otherwise a method that simulates individual states is to be preferred. 

Where the main interaction of the open system is with a measuring apparatus, which 
records the quantum state by counting, as with a photon counter, then the relative state or 
quantum jump method is the natural one and the best to use for numerical computations. 
The application of QSD produces rapid changes in state that look like quantum jumps 
p5| . p4| , but this is numerically inefficient. 

When measurement is not the main interaction with an environment, then QSD is often 
more efficient. Examples are the types of thermal interaction that are often represented 
by heat baths, interactions with a radiation field, and also continuous interactions with 
condensed matter or gaseous environments that are not to be considered as part of the 
system, as in line broadening, molecular dynamics, noise in quantum circuits, etc. In these 
cases it is often simplest to use QSD with a fixed basis, and this is likely to be more efficient 
numerically than jumps. 

Where the interaction with the environment is sufficiently strong to produce strong lo- 
calization during a significant part of the evolution, then MQSD is to be preferred. This is 
a very common situation, as the greater the environmental interaction, the greater the lo- 
calization becomes. MQSD is feasible where there are several degrees of freedom, and there 
would be no hope of using other simulation methods. In some limiting cases semiclassical 
methods can be used instead. MQSD is likely to be valuable in the those situations that lie 
between those that are sufficiently simple to be simulated by other numerical means, and 
those for which a semiclassical method can be used. 

Monte Carlo trajectory methods for open systems have been developed largely in con- 
nection with quantum optics, but they have been extended using QSD and MQSD to 
particles in traps []39[| , while Spiller and his collaborators have carried out some inter- 



esting QSD studies of quantum circuit elements with a view to applications to Josephson 
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junctions. 

However, the potential range of useful applications is far greater. In particular there has 
been no application to quantum effects in the dynamics of molecules which interact strongly 
with gases or with liquids. For isolated molecules, there are often transitions between many 
potential energy surfaces, and no clear way of localizing the molecule on one surface out of 
many. But in QSD, the interaction with the enviroment localizes the molecule into a small 
region around a phase point on one many possible surfaces, so that MQSD for a molecule 
interacting with its environment would have considerable advantages over the numerical 
solution of the Schrodinger equation for an isolated molecule. 

It is also clear that the noise in quantum circuits is an environmental effect that could 
well be simulated by MQSD. This applies even more strongly to quantum computation, 
where the potentially rapid decoherence produced by thermal interaction could be a crucial 
limitation. In these cases a lot could be learned from individual MQSD runs, without the 
need to make the large number of simulations required to get good statistics. It is not 
even clear that the density operator is of much relevance to such applications, for which 
the individual system is all- important. 
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Figure 1: Poincare surface for expectation values of x versus p for the forced, damped 
Duffing oscillator, plotted at times of constant phase t n = 2im. This is in the chaotic 
regime, with constants g = 0.3 and T = 0.125, and scaled up by a factor of 100 in x and p, 
effectively reducing h by a factor of 10 4 . The system is initially in the ground state. Note 
that it would take ~ 50, 000 states to represent this system with a non- moving basis. The 
cutoff probability is e = 5 x 10~ 3 . 



Figure 2: Photon number versus dimensionless scaled time r = K\t for a single second 
harmonic generation trajectory. Solid line: fundamental mode (ajai); broken line: second 
harmonic mode (ajj^). The parameters are k 2 / H\ = 0.25, 5i/ki = 5 2 / H\ = — 1, xl K i — 0-5, 
and f / Ki = 62. At time t — 0, the system is in the vacuum state. These parameters lie in 
the chaotic regime of the corresponding classical system. The cutoff probability is e = 10~ 3 . 
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